
// Store observed data in matrix for MATA
mata
mata clear

void mse_function_2dice(todo, theta, y, g, negH) {
    // Load data from Stata within function
    observed_data = st_data(., "reports")
    roll1_data = st_data(., "roll1")
    roll2_data = st_data(., "roll2")
    // Calculate one_eta values for first die
    one_eta = J(6, 1, .)
    one_eta[1] = 2/theta 
    for (i = 2; i <= 6; i++) {
        one_eta[i] = 2/((i-1)*(theta-2) + theta)
    }
    
    // Constants
    H_chance = 1/6
    
    // Calculate lying up fractions for first die (from highest to lowest)
    lu_frac = J(6, 1, .)
    lu_frac[6] = one_eta[6]
    for (i = 5; i >= 1; i--) {
        product = 1
        for (j = i+1; j <= 6; j++) {
            product = product * (1 - one_eta[j])
        }
        lu_frac[i] = one_eta[i] * product
    }
    
    // Calculate honest fractions for first die
    H_frac = J(6, 1, .)
    H_frac[6] = 1
    for (i = 5; i >= 1; i--) {
        H_frac[i] = 1 - sum(lu_frac[i+1::6])
    }
    
    // Calculate predicted probabilities for each observation
    predicted = J(rows(observed_data), 1, .)
    for (obs = 1; obs <= rows(observed_data); obs++) {
        r1 = roll1_data[obs] + 1  // Convert 0-5 to 1-6 for indexing
        r2 = roll2_data[obs] + 1  // Convert 0-5 to 1-6 for indexing
        
        // First die probability from DD model
        first_die_prob = H_chance * H_frac[r1] + lu_frac[r1] * (r1-1)/6
        
        // Second die is always uniform (1/6)
        predicted[obs] = first_die_prob * (1/6)
    }
    
    // MSE calculation
    residuals = (observed_data - predicted) * 100
    y = sum(residuals:^2) / (36)  // normalize
}

// Optimization
S1 = optimize_init()
optimize_init_technique(S1, "nr")
optimize_init_which(S1, "min")
optimize_init_evaluator(S1, &mse_function_2dice())
optimize_init_params(S1, 3)

theta_hat = optimize(S1)
mse_val = optimize_result_value(S1) 

st_numscalar("theta_hat", theta_hat)
st_numscalar("mse", mse_val)

theta_hat
mse_val
sqrt(mse_val)

end

display "Estimate of Theta is " theta_hat
display "corresponding MSE is " mse

// Calculate implied spread using estimated theta
scalar theta = theta_hat

// Calculate fractions for first die
forvalues i = 0/5 {
    local one_eta_`i' = 2/(theta + `i'*(theta-2))
}

// Lying up fractions for first die
local lu_frac5 = `one_eta_5'
local lu_frac4 = `one_eta_4' * (1-`one_eta_5')
local lu_frac3 = `one_eta_3' * (1-`one_eta_4') * (1-`one_eta_5')
local lu_frac2 = `one_eta_2' * (1-`one_eta_3') * (1-`one_eta_4') * (1-`one_eta_5')
local lu_frac1 = `one_eta_1' * (1-`one_eta_2') * (1-`one_eta_3') * (1-`one_eta_4') * (1-`one_eta_5')
local lu_frac0 = `one_eta_0' * (1-`one_eta_1') * (1-`one_eta_2') * (1-`one_eta_3') * (1-`one_eta_4') * (1-`one_eta_5')

// Honest fractions for first die
local H_frac5 = 1
local H_frac4 = 1 - `lu_frac5'
local H_frac3 = 1 - `lu_frac5' - `lu_frac4'
local H_frac2 = 1 - `lu_frac5' - `lu_frac4' - `lu_frac3'
local H_frac1 = 1 - `lu_frac5' - `lu_frac4' - `lu_frac3' - `lu_frac2'
local H_frac0 = 1 - `lu_frac5' - `lu_frac4' - `lu_frac3' - `lu_frac2' - `lu_frac1'

// Constants
local H_chance = 1/6

// Generate predicted values and display comparison
gen predicted = .
gen first_die_prob = .

forvalues i = 0/5 {
    local fdp = `H_chance' * `H_frac`i'' + `lu_frac`i'' * `i'/6
    replace first_die_prob = `fdp' if roll1 == `i'
}
replace predicted = first_die_prob * (1/6)

// Display first 12 combinations for verification
display "Implied vs Observed (first 12 combinations):"
list roll1 roll2 predicted reports in 1/12, clean noobs




rename predicted predicted_dd_$treat
drop  first_die_prob 

